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ABSTRACT 

We study the growth of massive black holes (BH) in galaxies using smoothed particle 
hydrodynamic simulations of major galaxy mergers with new implementations of BH 
accretion and feedback. The effect of BH accretion on gas in its host galaxy is modeled 
by depositing momentum at a rate ~ tL/c into the ambient gas, where L is the 
luminosity produced by accretion onto the BH and r is the wavelength-averaged optical 
depth of the galactic nucleus to the AGN's radiation (a free parameter of our model). 
The accretion rate onto the BH is relatively independent of our subgrid accretion 
model and is instead determined by the BH's dynamical impact on its host galaxy: 
BH accretion is thus self-regulated rather than "supply limited." We show that the final 
BH mass and total stellar mass formed during a merger are more robust predictions of 
the simulations than the time dependence of the star formation rate or BH accretion 
rate. In particular, the latter depend on the assumed interstellar medium physics, 
which determines when and where the gas fragments to form star clusters; this in turn 
affects the fuel available for further star formation and BH growth. Simulations over 
a factor of ^ 30 in galaxy mass are consistent with the observed Mbh — o relation 
for a mean optical depth of r ~ 25. This requires that most BH growth occur when 
the galactic nucleus is optically thick to far-infrared radiation, consistent with the 
hypothesized connection between ultra-luminous infrared galaxies and quasars. We 
find tentative evidence for a shallower Mbh — c relation in the lowest mass galaxies, 
a < lOOkms"^. Our results demonstrate that feedback-regulated BH growth and 
consistency with the observed Mbh — cr relation do not require that BH feedback 
terminate star formation in massive galaxies or unbind large quantities of cold gas. 

Key words: galaxies: evolution - galaxies: active 



1^ ■ 1 INTRODUCTION 



Feedback from an active galactic nucleus (AGN) has 
been invoked to resolve a number of observational prob- 
lems in galaxy formation: ( 1) to explain the tight ob- 
served i Ferrarese fc MerrittI |2000| : iGebhardt et"al1 l2000l : 
iHaring fc Rix 200^)" correlations between central black hole 
(BH) and galaxy properties such as the Mbh—ct and Mbh — 
M, relation s and the BH "fundamental plane" (ISilk fc ReesI 
19981: iKindlioO S; Murr ay et al.ll2005l:lDi Matteo et al.ll2005l : 



Sazonov et 811120051 : iHopkins et al.| |2007l '). (2) to shut off star 
formation in elliptical galaxies (e.g., by blowing gas out of 
the galajcy), ther eby explaining how ell i pticals become "re d 
and dead" fe.g.. ISpringel et al.l l2005al : ICiotti et al.ll2010D . 
(3) to heat the hot intracluster plasma (ICM) in groups 
and clusters, thereby suppressing coolin g and star f orma- 
tion in these environments (e.g.. Tabor fc BinnevI 1 19931 : 
ICiotti fc Ostrikeill 19971 : ICroton et al.ll2006h . and (4) to help 
explain "cosmic downsizing," namely the fact that both star 



formation and AGN activity resi de in progressively lowe r 
mass halos at lower redshifts (e.g.. lScannapieco et al.|[2005l ). 

It is plausible that AGN perform the roles desired 
of them, but this is by no means certain. Understanding 
whether this is indeed the case requires developing more so- 
phisticated theoretical models that can be compared quan- 
titatively to observations. There are several key theoretical 
problems that must be addressed in order to better under- 
stand the role of massive BHs in galaxy formation, and to 
understand the properties of massive BHs themselves. The 
first is the problem of AGN fueling, i.e., how is gas trans- 
ferred from galactic scales (~ 0.1 — 1 kpc) to the vicinity 
of the massive BH (< 0.1 pc)? A second key problem is the 
problem of AGN feedback: how do energy and momentum 
generated by accretion onto a central BH - in the form of 
radiation and outflows - couple to the surrounding gas, and 
how does this affect star formation and the growth of the 
BH itself? 
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Much of the recent work addressing the impact of BHs 
on g alaxy formation has u s ed quahtatively simila r physics 
(e.g., ISprineel et al]l2005bl : I Johansson et al.ll2009l ). For ex- 
ample, many calculations assume that a BH of mass Mbh 
will accrete mass at a rate proportional to the Bondi rate 
l|Bondi, |1952): 



(1) 



where p is the density of the surrounding g the 
sound speed of that gas, and / ~ 10 — 100 is a factor 
taking into account the possible multi-phase structure of 
the gas and tha t the sphere of influen ce of the BH is of- 
ten not resolved ([Booth fc Schav3l2009h . There is, however, 
little justification for using equation [T] The Bondi accre- 
tion rate estimate assumes that the gas surrounding the BH 
is spherically symmetric. When the gas is not spherically 
distributed, the rate of angular moment um transport de ter- 
mines the BH accretion rate (e.g.. ,Shlosman et al.lll990l V It 
is generally believed that the progenitors of todays > L* el- 
lipticals are gas-rich disk galaxies, the mergers of which lead 
to lu minous starbursts and the growth of the ce ntral massive 
BHs (jSanders etallllQSSl : [Hopkins et all bOOSV Most of the 
gas in disk galaxies, merging g alaxies, lumino us starbursts 
([Pownes fc Solomon|[l9 98: Tacc oni et al.|[200^ '). and nearby 
luminous AGN ([Ho et al.ii2008l ') appears to reside in a rota- 
tionally supported disk. There is thus no reason to expect 
that the spherically symmetric Bondi rate provides a good 
estimate of the BH accretion rate in gas rich galaxies. Even 
in the central ~ parsec of own galaxy, where the ambient 
gas is hot and pressure supported, the Bondi accretion rate 
fails by orders of m agnitude to p r edict the accretion rate 
onto the central BH ([Sharma et aLi l2007[ ). 

There are a number of way s that an AGN can s trongly 
influence its surroundings (e.g.. [Ostriker et al.ll201(il ^. Rela- 
tivistic jets inject energy into intracluster plasma and may 
be the primary mechanism suppressing co oling flows in 
galaxy clusters ([McNamara fc NulsenI [20071 ). even though 
the details of how the energy in the jet couples to the 
plasma in a volume filling way are not fully understood 
llVernaleo fc Revnoldsl20"06f ). On galactic scales, a wind from 
an accretion dis k around t he BH can drive gas out of 
the galaxy (e.g., [Kind \200^ ) as could cosmic-ray proton s 
produced by a radio loud AGN ([Sironi fc Socrates| [2010t ). 
In addition, the AGN's radiation can strongly affect the 
surrounding gas, bo th by Compton heating/cooling (e.g., 
ISazonov et al.l |2005| ) and by the mome ntum imparted as 
UV radiation is absorbed by dust grains l[Chang et al.ll 19871 : 
[Sanders et al] [l988[ : [Murray et aT|[2005[ ) . 

This diversity of feedback mechanisms can be roughly 
separated into two broad classes: energy and momentum 
injection. We believe that momentum injection is the domi- 
nant mode of feedback for most of the gas in a galaxy, largely 
because of the very short cooling times of dense gas. For ex- 
ample, if a BH radiates at ^ 10*® erg s~^ with a typical 
quasar spectrum, only gas with n < 1 cm~"^ can be heated 
to the Compton temperature within ~ 100 pc. However, the 
mean gas densities in the central ~ 0.1 — 1 kpc of luminous 
star forming galaxies are ~ 10^~^cm~^ ([Pownes fc SolomonI 
1 19981 : iTacconi et al.l [2006l 'l. At these densities, the cooling 
time of gas is sufficiently short that it is unable to retain 
much injected energy - be it from the AGN's radiation or 



from shocks powered by AGN outflows. Thus it is largely the 
momentum imparted by AGN outflows and by the absorp- 
tion and scattering of the AGN's radiation that dominates 
the impact of the AGN on dense gas in galaxies. Since it 
is the dense gas that fuels star formation and the growth 
of the BH itself, it is critical to understand the impact of 
momentum feedback on this gas|3 

In this paper, we present simulations of major mergers 
of spiral galaxies using a model for the growth of BHs that 
includes (1) a BH accretion rate prescription motivated by 
the physics of angular momentum transport and (2) AGN 
feedback via momentum injection (e.g., radiation pressure). 
Some results of this model appear in a companion Letter 
([PeBuhr et al.[[2009[ ). The remainder of this paper is orga- 
nized as follows. Section [5| presents a summary of our meth- 
ods, including a description of the model galaxies ( i]2.ip . the 
model for star formation and the interstellar medium ( i]2.2|) . 
our BH accretion and feedback model (' >i2.3|l and a summary 
of our parameter choices (Sj2]4}. Section [3| shows the results 
of applying this model to BH growth and star formation in 
major mergers of gas-rich galaxies. In section[3|we show that 
our model of BH growth and feedback produces a reasonably 
tight Mbh — o" correlation similar to that observed. Finally, 
in section[S|we discuss our results and compare our approach 
to previous models in the literature. Appendix [X| presents 
resolution tests for our fiducial simulation while Appendix 
[B] presents some of the tests used to verify the BH accretion 
and feedback models that we have implemented. 

2 METHODOLOGY 

W e use a non-pub Uc update of the TreeSPH code GAPGET- 
2 l[Springel[[2005 ') provided by V. Springel to perform sim- 
ulations of equal-mass mergers of galaxies. This version 
of the code includes t he eff ective star formation model of 
[Springel fc HernquistI ( [2003[ ) but contains no AGN feedback 
physics. We modified the code further to implement models 
for massive BH growth and AGN feedback. The details of 
the simulations are described in the following subsections. 
The Appendices present resolution tests and some of the 
tests we performed to verify our implementation of the BH 
accretion and feedback model. 

2.1 Initial Conditions and Galaxy Parameters 

Each model galaxy u sed in our ma j or mer ger simulations 
is similar to those in [Springel et all ([2005bl ). They include 
a spherical halo of coUisionless dark matter, a centrifugally 
supported disk of gas and stars, a stellar bulge, and a central 
point mass representing a black hole. The code used to gen- 
erate the initial conditions was provided by V. S pringel and 
is identical to that used in lSpringel et al.l l[2005bf ) except for 
one change that will be described below. 

Table 1 lists the relevant galaxy and simulation parame- 
ters for the key merger simulations we focus on in this paper. 
The simulations are all major mergers of equal mass galax- 
ies. The fiducial simulation (top entry) assumes a mass of 
1.94 X IO^^Mq for each merging galaxy, of which 4.1% is 



^ These conclusions do not apply to dilute plasma in the intra- 
cluster or intragroup medium. The densities there are sufficiently 
low that the plasma can be efficiently heated by an AGN. 
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assigned to the gas and stars in the disk, 1.36% is assigned 
to the stars in the bulge, and the rest is in a dark matter 
halo. The initial mass fraction of gas in the disk is fg = 0.1. 
This run uses a total of A'p = 1.6 x 10^ particles with 6 x 10'' 
dark matter particles, 2 x 10^ particles each in the gaseous 
and stellar disk, and 10^ particles for the stellar bulge. This 
run has a Plummer equivalent gravitational force softening 
of e = 47 pc. 

To test the dependence of the results of our fiducial 
simulation on the model and simulation parameters, we have 
run a number of additional simulations, varying the gas frac- 
tion (fg = 0.3 vs 0.1), bulge-to-disk mass ratio (0.2 vs 0.33), 
total galaxy mass (from 0.1 to 3 of the fiducial value), sim- 
ulation particle number (from A'p = 1.6 x 10^ to 2.4 x lO*"), 
force softening (e — 22 to 102 pc), as well as the parameters 
in the black hole model ( descr ibed in § 2.4 below). 

We use a lHernguistl (|l990l ) density profile for the struc- 
ture of the dark matter halo: 



Phalo{r) 



Mhalo a 

2n r{r -f a)^ 



(2) 



The scale length a of the halo is set by requiring that the 
halo enclose the same mass within the virial radius as an 
NFW profile, and that the densities match at small radii. 
These conditions yield a relationship among the halo scale 
length, a, the corresponding NFW s cale length, rs, and 
the concentration of the NFW halo, c l|Navarro et al.|[l996l : 
ISpringel et al. l l2005bD : a = r42[ln [l + c)- c/(l + c)\Y'^ . 
The halos used in this work all have a concentration of c = 9. 

The stellar and gaseous disks both initially have expo- 
nential surface density profiles: 



E(i?) 



2'kRI 



exp - — 



(3) 



where Mi is the total mass of the component of interest and 
Rd is the disk scale length, which is initially the same for the 
stellar and gaseous disks. The disk scale length for the fidu- 
cial simulation is Rd = 3.5 kpc, which corresponds to the 
disk having approximately the same angular momentum per 
unit mass as a halo with a spin parameter of 0.033. For simu- 

1/3 

lations with different disk masses, we use Rd oc which 
is consistent with the observed relation (Shc n et al.ll200v? l. 
The stellar disk's vertical structure is given by the standard 
sech^ {z/zo) profile, where the vertical scale height zq is ini- 
tially set to zo = -Rd/5 at all radii. Unlike the stellar disk, 
the gaseous disk's vertical structure is determined by hydro- 
static equilibrium given the assumed sound speed/equation 
of state of the gas (discussed below). Setting up this initial 
vertical hydrostatic equil ibrium require s an iter ative proce- 
dure that is described in ISoringel et al.l ()2005bl l . 

The stellar bulges also have Hernquist density profiles. 
The scale length of the bulge Rb is specified as a fraction of 
the disk scale length, Rd- In the fiducial simulation, Rt, = 
Rd/5. For different bulge masses, we use the scaling relation 

1/2 

Rb oc , which is motivat ed by the observed mass-radius 
relation of elliptical galaxies (|Shen et al.ll200j l. 

In our simulations, two galaxies with identical struc- 
ture are placed on a prograde orbit. For simulations at our 
fiducial mass of 1.94 x W^^Mq (for each galaxy), the initial 
separation of the two galaxies' centers is 142.8 kpc. The or- 
bit has approximately zero total energy, which corresponds 
to an initial velocity for each galaxy of 160 km s~^; the ve- 



locity is directed at an angle of 28 degrees from the line 
connecting the centers of the two galaxies. In order to break 
the symmetry of the problem, the individual spin axes of 
the galaxies have a relative angle of about 41 degrees, with 
one galaxy of the pair having an inclination with respect to 
the orbital plane of 10 degrees. For the simulations with dif- 
ferent overall masses, the orbital parameters are scaled by 
M^^^ , so that the time to first passage and the time to final 
merger are similar to those in the fiducial run. 

2.2 Interstellar Medium Model 

The version of GADGET we use includes 
ISpringel fc Hernguistl {2003) 's sub-resolution model for 
the interstellar medium (ISM). This model treats the gas as 
a two phase medium of cold star forming clouds and a hot 
ISM. When cooling and star formation are rapid compared 
to the timescale for adiabatic heating and/or cooling (which 
is nearly always the case in our calculations), the sound 
speed of the gas is not determined by its true temperature, 
but rather by an effective sound speed that averages over 
the multi-phase ISM, turbulence, etc. The effective sound 
speed as a function of density can be interpolated freely 
between two extremes using a parameter qaos- At one 
extreme, the gas has an effective sound speed of lOkms"^, 
motivated by, e.g., the observed turbulent velocity in 
atomic gas in nearby spirals; this is the "no-feedback" 
case with Qeos = 0. The opposite extreme, Qcos ~ 1, 
represents the "maximal fee dback" sub-resolution model of 
ISprineel fc Hernauisd tOO±. mo t ivated by the multiphase 
ISM model of iMcKee fc Ostrikeil (|l977l ): in this case, 100% 
of the energy from supernovae is assumed to stir up the 
ISM. This equation of state is substantially stiffer, with 
effective sound speeds as high as ~ 200kms~^. Varying 
geos between these two extremes amounts to varying the 
effective sound speed of the ISM, with the interpolation 



Cs = V (?cos ci [q = 1] + (1 - gcos) ci [q = 0] 



(4) 



In addition to this effective equation of state, GAD- 
GET models star formation by stochastically converting gas 
particles into star particles at a rate determined by the gas 
density. 



PSF = 



1-P 1/2 



oc p 



3/2 



(5) 



where /3 = 0.1 is the fraction of the mass of a stellar popula- 
tion returned to the ISM by stellar evolution. The parameter 
is the characteristic timescale for gas to be converted into 
stars at the threshold density pth = 0.092 cm~^; pc ~ p 
is the density of the cold clouds, which is related to the 
density of the SPH particle by equations (17) and (18) of 
ISpringel fc Hernauisd (|2003l ). For a given gas equation of 
state, the parameters in equation [5] can be adjusted to pro- 
duce a global star formation la w similar to the observed 
Kennicutt-Schmidt relations (Sp ringel et al.|[2005bl ). 

For parameters in the equati on of state model th at 
have been used in previous work (ISpringel et al.ll2005bD - 
TsN = 4 x 10** K, ^0 = 4000, tl = 8.4 Gyr and qsos = 0.5 - 
we find that the model overpredicts the sound speed relative 
to the observed "turbulent" velocities of galaxies, i.e., the 
non-t hermal line widths (see Fig. 1 of iHopkins fc Quataerd 
I2OO9I for a compilation of relevant data). For instance, the 
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Table 1. Simulation Parameters 



Run Name 


Mtot 


fafi 




Np 


e 


Race 
€ 












[106] 


[pc] 




fid 


1.0 


0.1 


0, 


.33 


1.6 


47 


4 


fidNof' 


1.0 


0.1 


0, 


.33 


1.6 


47 


4 


fidSa 


1.0 


0.1 


0, 


.33 


1.6 


47 


4 


fid6a 


1.0 


0.1 


0, 


.33 


1.6 


47 


4 


fidTau 


1.0 


0.1 


0, 


.33 


1.6 


47 


4 


fidt25 


1.0 


0.1 


0, 


,33 


1.6 


47 


4 


fidSeps 


1.0 


0.1 


0, 


,33 


1.6 


47 


8 


fidafg 


1.0 


0.1 


0, 


,33 


1.6 


47 


4 


fidq2'* 


1.0 


0.1 


0, 


,33 


1.6 


47 


4 


fidqOT'^ 


1.0 


0.1 


0, 


,33 


1.6 


47 


4 


big 


3.0 


0.1 


0, 


,33 


1.6 


68 


4 


big6a 


3.0 


0.1 


0, 


,33 


1.6 


68 


4 


mid 


0.3 


0.1 


0, 


,33 


1.6 


32 


4 


small 


0.1 


0.1 


0, 


,33 


1.6 


22 


4 


small6a 


0.1 


0.1 


0, 


,33 


1.6 


22 


4 


smaUqO?*^ 


0.1 


0.1 


0, 


,33 


1.6 


22 


4 


fg 


1.0 


0.3 


0, 


,33 


2.4 


47 


4 


smallfg 


0.1 


0.3 


0, 


,33 


2.4 


22 


4 


bulge 


1.0 


0.1 


0, 


,20 


1.6 


47 


4 


LRfid 


1.0 


0.1 


0, 


,33 


0.16 


102 


4 


MRfid 


1.0 


0.1 


0, 


,33 


0.48 


70 


4 


MRfidNof'' 


1.0 


0.1 


0, 


,33 


0.48 


70 


4 


LRfidNof'' 


1.0 


0.1 


0, 


,33 


0.16 


102 


4 


fidvol 


1.0 


0.1 


0, 


,33 


1.6 


47 


8.62 


MRfidvol 


1.0 


0.1 


0, 


,33 


0.48 


70 


5.97 



a 


r 


, new 




Mbh,p 


(Tf 






y^\j i J J 


[lO^Mnl 

1^ J. vy ^ V J. J 




[kms^i] 


0.05 


10 


1.34 


1.49 


1.33 


169 


0.15 





1.36 


18.1 


13.5 


170 


0.15 


10 


1.34 


1.03 


0.90 


168 


0.3 


10 


1.35 


0.86 


0.77 


167 


0.05 


3 


1.36 


5.05 


4.31 


163 


0.05 


25 


1.35 


0.39 


0.35 


169 


0.05 


10 


1.35 


2.70 


1.76 


163 




10 


1.32 


1.21 


1.02 


169 


0.05 


10 


1.30 


1.40 


1.16 


168 


0.05 


10 


1.32 


1.52 


1.36 


164 


0.05 


10 


3.08 


6.24 


5.27 


232 


0.3 


10 


4.17 


7.86 


5.15 


227 


0.05 


10 


0.39 


0.38 


0.26 


115 


0.05 


10 


0.13 


0.24 


0.13 


82 


0.3 


10 


0.13 


0.25 


0.24 


84 


0.05 


10 


0.12 


0.06 


0.05 


81 


0.05 


10 


4.41 


7.10 


5.53 


159 


0.05 


10 


0.36 


0.31 


0.23 


98 


0.05 


10 


1.38 


1.44 


1.25 


161 


0.05 


10 


1.34 


1.65 


0.93 


164 


0.05 


10 


1.35 


2.92 


2.40 


168 


0.15 





1.34 


13.5 


11.4 


167 


0.15 





1.31 


13.1 


11.4 


175 


0.05 


10 


1.39 


3.22 


2.45 


164 


0.05 


10 


1.36 


3.30 


1.92 


164 



Columns are defined as follows: Mtot is the total mass in the simulation, fg^a is the initial gas fraction of the disk, Mf^/M^ is the bulge 
to disk mass ratio, A^p is the total number of particles used in the simulation, e is the Plummer equivalent gravitational force softening, 
cc^ and T are the parameters of the BH accretion and feedback model f ^2.311 . M^^new is the total mass of new stars formed during the 
simulation, Mbhj and Mbh,p S't'e the masses of the BH at the end of the simulation and after the peak of accretion (defined to be when 
the accretion rate drops to one tenth its maximum value), respectively, and o"/ is the stellar velocity dispersion of the merger remnant (Q. 



a Mfii = 3.88 X IO^^Mq. 

These runs had no AGN feedback. 

a was set by the gas fraction within Race using a = 3/^ . 

I 

" The ISM equation of state is defined using q^os = 0.2 (see N2.2I I. 

^ The ISM equation of state is defined using Qcos — 0.07 (see [[g^J. 

above model parameters imply Cs ~ 30 km at n ~ 1 
cm~^ and Cs ~ 110 km s^^ at n ~ 10^ cm~^. These val- 
ues are too large by a factor of ~ 2 — 3 compared to the 
random velocitie s inferred from atomic and molecular line 
observations (Downes fc Solomon . 1998). To account for this, 
we set Qeos = 0.5 and then modified GADGET by reducing 
the pressure everywhere by a factor of 10. This reduces the 
effective sound speed by a factor of ~ 3 and is thus more 
consistent with observations. This reduction in ISM pres- 
sure is also used in the initial conditions when setting up 
vertical hydrostatic equilibrium for the gas. Changing the 
pressure requires changing the equation of state parameters 
to TsN = 6.6 X 10** K, Ao = 6600, and t° = 13.86 Gyr to 
maintain an average star formation rate of 1 Mq yr~^ for 
an isolated galaxy with our fiducial Milky Way like mass. 
In ^3.31 we compare our fiducial calculations with this re- 
duction in pressure to models with smaller values of Qcos, 
0.07 and 0.2; these also have smaller "sound speeds" more 
comparable to the observed random velocities of galaxies. 



The reduction in the sound speed decreases the Jeans 
length and mass, making it numerically more prohibitive to 
resolve these critical scales. For the simulations presented 
here, we are careful to use sufficient numbers of particles so 
that the Jeans length and mass are always adequately re- 
solved. The higher gas fraction simulations require a higher 
particle number as a result (see Table [TJ. The reduction in 
sound speed also makes it more likely that the gas will frag- 
ment by gravitational instability into clumps (ala molecular 
clouds), as we shall discuss in detail later. This fragmenta- 
tion is real, not numerical; artificially increasing the sound 
speed to eliminate it is not necessarily physical and could 
give incorrect results. On the other hand, we do not include 
sufficient physics in our ISM model to describe the forma- 
tion and disruption of molecular clouds so our treatment of 
the resulting clumping is also not correct. In i|3.3l we discuss 
which of our results are the most sensitive to uncertainties 
related to local gravitational instability in the ISM. 
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2.3 Black Hole Accretion and Feedback 

2.3.1 Black Hole Accretion Model 

We include a BH as an additional collisionless particle at 
the center of each galaxy. We model the accretion of the 
surrounding gas onto the BH, via the transport of angular 
momentum, using 

Mvisc — SvraE-^ (6) 

where S is the mean gas surface density, Q, is the rotational 
angular frequency, and a is the dimensionless viscosity (a 
free parameter of our model). We compute E and Cs by tak- 
ing an average of the properties of the SPH particles in a 
sphere of radius Race centred on the BH. The radius Race is 
typically set equal to four times the gravitational force soft- 
ening length, i.e.. Race = 4e, although we explore alternate 
choices as well. We find that estimating the rotation rate 
using ~ GM{< Race) /Race is more numerically robust 
than actually calculating the rotation and angular momen- 
tum of the gas particles within Race. 

Althou gh equation (El) is reminisc ent of the alpha pre- 
scription of IShakura fc SunvaevI (119731 ) . in our formulation 
a characterizes not only the efficiency of angular momen- 
tum transport, but also the uncertainty due to the fraction 
of the inflowing gas that is turning into stars vs. being ac- 
creted onto the central BH. The physical mechanisms driv- 
ing gas from ~ kpc to ~ 0.1 pc are not fully understood, but 
non -axisymmetric gravitational torques are likely re sponsi- 
ble l|Shlosman et al.lll989l : lHopkins fc Quataertll2009l ). Using 
numerical simulations t hat focus on the nuc le i of g alaxies 
(from ~ 0.1 - 100 pc) iHopkins fc QuataertI (|2009l ) simu- 
late the conditions under which there is significant gas in- 
flow to < 0.1 pc. They argue that the net accretion rate is 
not a strong function of the gas sound speed (unlike both 
eqns [T] and [5| because non-axisymmetric gravitational per- 
turbations produce orbit crossing and strong shocks in the 
gas. The resulting inflow rate depends primarily on the non- 
axisymmetry in the potential, rather than the thermody- 
namics of the gas. Nonetheless, equation (|S| evaluated at 
~ 100 pc and with a 0.1 approximates the accretion rate 
at small radii in their simulations, albeit with substantial 
scatter (factor of ~ 10). Given that one of our key results 
discussed in ^is that the accretion rate is not sensitive to 
the exact value of a, we believe that equation ^ is sufficient 
for the exploratory calculations in this paper. 

2.3.2 Mass of the Black Hole Particle 

In our galaxy merger simulations, the two BHs are initially 
far apart but approach each other in the late stages of the 
merger. When the BH particles have a separation of less than 
Race we consider them to have merged. When this occurs, 
we sum the individual masses of the two BH particles and 
set one of the particles to have this mass. This particle is 
then moved to the center of mass of the two BH system and 
given the velocity of the center of mass frame. The other BH 
particle is removed from the region. 

The BH particles are subject to stochastic motion due 
to interaction with the stellar and gaseous particles, which 
leads to inaccuracy in the position of the BH and noise in 
the estimate of the accretion rate. To reduce this numeri- 
cal "Brownian" motion, the BH particles are given a large 



20 

scaled with the overall mass for other simulations. As a re- 
sult, the BH particle is a factor ~ 100 more massive than 
the halo particles, and a factor ~ 10'' more massive than the 
stellar and gaseous particles. We artificially increase the BH 
particle mass solely to reduce numerical relaxation effects. 
This does not result in spurious dynamical effects on the 
central stars, gas, and dark matter since the BH's sphere of 
gravitational influence extends to < 10 pc for the flducial 
simulation, which is signiflcantly smaller than our typical 
force softening of ~ 50 pc. 

For the results presented below, the "real" mass of the 
BH (= Mbh) is computed by integrating the accretion rate 
of equation (|6]) in time. The gas particles are not removed 
as the BH mass increases. Instead, the gas particles have an 
additional label that tracks whether or not they have been 
"consumed." We track how much mass the BH should have 
consumed via accretion at a given time, and the mass of 
gas that has been consumed. When there is a mis-match, 
we tag a number of gas particles within Race (chosen at 
random) as "consumed" until the total mass accreted by the 
BH is correct. Particles that have been consumed no longer 
contribute to the accretion rate estimate, even if they are 
inside Race. This implementation prevents any gas particle 
from providing more than its mass to the integrated mass 
of the BH. 



2.3.3 Feedback from the Black Hole 

In our simulations, the AGN is assumed to couple to the sur- 
rounding gas by depositing momentum into the gas, directed 
radially away from the BH. This crudely approximates 
the effects of (1) strong out flows and/or cosmic-ray pres- 
sure produced by the AGN l|Kinel l2003l : ISironi fc Socrates! 
2010) and (2) radiation pressure produced by the absorp- 
tion and scattering of th e AGN's radiation by dust in the 
ISM (|Murrav et al.ll2005l '). We focus on the latter when mo- 
tivating the parameters used in our models. 

To accurately account for the impact of the AGN's ra- 
diation on gas in its host galaxy would require a radiative 
transport calculation, which is beyond the scope of the cur- 
rent work. Instead, we model this radiation pressure by de- 
positing a total momentum per unit time of 



p = where L = min (rjMyiseC^ , Lsdd^ (7) 

radially away from the BH into the SPH particles within 
a distance of Race of the BH particle. This momentum is 
equally distributed among the particles so that each par- 
ticle experiences the same acceleration. We use a radiative 
efficiency of = 0.1 in all simulations. The physical picture 
behind our feedback model in equation ((Tjl is that the feed- 
back is produced by the absorption of the ultraviolet light 
from the AGN by dust in the surrounding gas, and the sub- 
sequent reemission of infrared radiation that must diffuse 
its way out of the nuclear region. As described shortly, the 
parameter r is the total infrared optical depth of the nuclear 
region. 

To motivate equation (O in more detail, we note that 
AGN radiate most of their radiation in the ultraviolet. The 
opacity of dusty gas to UV radiation is kuv ^ 10^ cm^ g~^, 
so that only a surface density of --^ 10~^ g cm~^ is required 
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to absorb the UV radiation. This is far less than the typ- 
ical radial column density of gas in the central ~ 0.1 — 1 
kpc of luminous star forming galaxies, galaxy mergers, or 
our simulations (see Fig. [2] below). As a result, the UV ra- 
diation is efficiently absorbed, except perhaps along polar 
lines of sight. The absorption and scattering of the UV ra- 
diation deposits a momentum per unit time of L/c into the 
ambient gas, assuming for simplicity that all of the UV ra- 
diation is absorbed. If the infrared optical depth is > 1, the 
infrared radiation re-emitted by the dusty gas must diffuse 
out through the nuclear region; doing so deposits an addi- 
tional momentum per unit time of tL/c, where r ~ hirYI is 
the infrared optical depth and kjr ~ few- 10 cm^ g~^ is the 
infrared opacity for the radiation temperatures of interest 
~ 100 - 1000 K. The net force due to the UV and infrared 
radiation is thus p ~ (1 -I- t)L/c ~ tL/c, i.e. equation 
for r > 1, which is valid in our calculations near the peak of 
activity when the BH gains most of its mass. 

In our calculations we use a constant value of r rather 
than a time variable r given by r = kirTi. Given the sim- 
plicity of our feedback model relative to a true radiative 
transfer calculation, this is not an unreasonable approxima- 
tion. It is also easier to isolate the effects of varying r when 
it is constant in time. 

As noted above, we apply the force in equation ((7} to 
all particles within a distance Race of the BH. A more ac- 
curate treatment would be to apply the force out to the 
point where the column is ~ «;7b> i-^-i to where the optical 
depth to infinity is ~ 1. At many times, however, this ra- 
dius is unresolved. Moreover, it is possible that the photons 
diffuse primarily along the rotation axis of the gas, rather 
than in the orbital plane. As a result, the radiation pressure 
force will be applied primarily at small radii. This is why we 
apply the force only within Race- One consequence of this 
is that the number of SPH particles experiencing the feed- 
back, A'^, will change as gas moves in and out of Race- Thus, 
the strength of feedback felt by an individual particle will 
change with time. However, because the SPH particles are 
coUisional, they readily share this momentum with neigh- 
boring gas particles. In test problems described in Appendix 
B the effects of our feedback model are essentially indepen- 
dent of N and Race- The results are not quite so clean in our 
full simulations (see i]3.2l and Appendix|3, but nonetheless 
none of our major results depend sensitively on the region 
over which the feedback force is applied. 

One might worry that if the number of particles within 
Race were too small, the momentum supplied to a single 
particle would become large enough to artificially acceler- 
ate the particle to the escape velocity. The minimum A'' re- 
quired to avoid this is actually quite modest for the range 
of luminosities in our calculations, and for the simulations 
presented here this concern is never an issue (although it is 
for some of the test problems in Appendix B). 

2.4 Parameter Choices for the Black Hole Model 

Our model for BH growth and feedback contains three free 
parameters: (1) a determines the magnitude of the accretion 
rate onto the BH; (2) r determines the total radiation pres- 
sure force produced by accretion onto the BH; it is roughly 
the optical depth to the far IR in the nuclear region; and (3) 
Race is the radius of the spherical region within which the 
accretion rate is determined and the feedback is applied. Our 



fiducial values for these parameters are a — 0.05, r = 10, 
and Race — 4e (where e is the gravitational force softening). 
We now motivate these particular choices. 

The fiducial value of the viscosity used in this work 
is a = 0.05, motivat ed by t he rough consistency between 
the resulting M and iHopkins fc Quataert (2009) 's numer- 
ical simulations of gas inflow from ~ 100 pc to ~ 0.1 pc 
(although there is fa ctor of ~ 10 scatter in the latter that 
is not captured here) . [Hopkins fc Quataerd (|2009l )'s calcula- 
tions in fact require a more complicated subgrid accretion 
model that depends on additional parameters such as the 
bulge to disk ratio of the galaxy (because this influences the 
strength of non-axisymmetric torques); this will be explored 
in more detail in future work. In addition to a = 0.05, we 
also carried out simulations with a = 0.15 and a = 0.3, and 
found no signiflcant differences, for reasons explained below. 

We use a constant value (with time) of r = 10 in most of 
our simulations. This is motivated by far infrared opacities 
of Km ^3 — 10 cm^ g~^ and surface densities of E ~ 1 — 10 g 
cm~^ within Race during the peak of activity in our simula- 
tions. These surface densities are also consistent with those 
directly measured in the nuclei o f ultra-luminous infrared 
galaxies l|Downes fc Solomon|[l9 9§). Given the uncertainties 
associated with the radiative transfer of far infrared photons 
in galactic nuclei, it is not possible to more accurately esti- 
mate the effective value of r without detailed radiative trans- 
fer calculations. As we shall demonstrate explicitly, however, 
the exact value of r is also not that critical for the qualita- 
tive effects of AGN feedback; the value of r does, however, 
strongly affect the final value of the BH mass. 

In choosing a value for Race, we must satisfy Race > e 
in order to avoid numerical artifacts. In addition, we find 
that the BH particle remains within 4e of the centre of mass 
of the system at nearly all times, but it can wander around 
within this region. As a result, 4e is the smallest we can 
make Race without having noise induced by the BHs motion. 
This choice corresponds to several hundred pc in our typical 
simulation. Larger values of Race are unphysical because (1) 
the accretion rate should only depend on the gas close to the 
BH; i.e., the transport of gas from, for example, ~ 8e to ~ 2e 
is presumably adequately described by our simulations so we 
should not try to also account for this in our subgrid model, 
and (2) the radiation pressure force produced by the AGN 
(and the re-radiated infrared photons) is likely concentrated 
at relatively small radii, for the reasons described in H2.3I 



3 GALAXY MERGER SIMULATIONS 

Table [T] summarizes the simulations we focus on in this pa- 
per, including the resolution, the parameters that specify 
the initial conditions for the merging galaxies, the param- 
eters that specify the BH accretion and feedback models, 
and the final properties of the merger remnants (stellar and 
BH mass and velocity dispersion). We begin by describing 
the results from our fiducial simulation (top row in Table [TJ 
and then discuss simulations that vary a single parameter 
of the feedback model relative to the fiducial run. We have 
also performed simulations at different overall galactic mass 
scales, initial gas fractions, and numerical resolution. The 
latter resolution tests are presented in Appendix A. 
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Figure 1. Top: The separation of the black hole particles as a 
function of time in the fiducial simulation. The blue circles la- 
bel the times of the images shown in Figure [3] Middle: The star 
formation rate as a function of time for the fiducial simulation 
(black) and for the run with no feedback (red; run fidNof). Bot- 
tom: The viscous accretion rate, M^isc (black), and Eddington 
rate (grey), as functions of time for the fiducial simulation. The 
critical Mc at which radiation pressure balances gravity (eq. [Sjl is 
shown within a radius of Race (red; solid). The increase in star 
formation and BH accretion after first passage {t ~ 0.75 Gyr) 
is due to the fragmentation and inspiral of large gaseous/stellar 
clumps (Fig. [Sj, while the much larger increase at final coales- 
cence is due to inflow of diffuse gas caused by non-axisymmetric 
torques. The latter physics dominates the total stellar and BH 
mass formed during the merger. 

3.1 The Fiducial Simulation 

The top panel of Fig. [T]shows the separation of the BH parti- 
cles for the fiducial simulation, while the middle panel shows 
the total star formation rate (in both galaxies) for simula- 
tions with (black) and without (red) BH feedback. The first 
close passage of the two galaxies is around t — 0.33 Gyr and 
the system then undergoes a few short oscillations as the 
BHs finally settle into a merged state around t — 1.65 Gyr. 
The star formation rate increases following the first pas- 
sage, with a much larger increase in the star formation rate 
during the final merger of the galaxies. The bottom panel 
of Fig. [T] shows the BH accretion rate determined from 
equation [6] (black) and the Eddington accretion rate (grey; 



Medd = Ledd/O-lc^); the initial BH mass is 1.4 x W' Mq but 
as long as it is not too large > 10* Mq, the precise initial BH 
mass is unimportant for our conclusions. In this and similar 
plots throughout the paper, the value of M plotted before 
the BHs merge is for the BH in the galaxy with the smaller 
initial inclination relative to the orbital plane; the BH accre- 
tion rate for the other galaxy is comparable to that shown 
here. The evolution of the accretion rate is similar in many 
of the simulations we have carried out, with an initial period 
of activity after the first passage of the merging galaxies, and 
another period of even higher M after the final coalescence 
of the galaxies and BHs. The latter active episode is when 
the merged BH gains most of its mass. In particular, the BH 
reaches the Eddington limit, allowing the mass of the BH to 
grow exponentially for a few hundred Myr. 

iDeBuhr et al] (|2009l ) showed that the BH accretion 
and feedback model presented in this work leads to self- 
regulated BH growth, due to a competition between the 
(inward) gravitational force produced by the galaxy as a 
whole and the (outward) r adiation pressure fo rce produced 
by the central AGN (eq.[7| (jMurrav et al.l2005h . For a spher- 
ically symmetric system, equating these two forces leads to 
tL/c = 4:fga*/G, where cr^ = GMt/2Racc, Mt is the to- 
tal mass inside Race, and we have evaluated these expres- 
sions within Race, where our accretion rate is determined 
and feedback is implemented. Equivalently, there is a criti- 
cal accretion rate Mc, analogous to the Eddington rate, at 
which the two forces balance: 

The bottom panel of Fig. [T]shows Mc for our fiducial simula- 
tion, evaluated within Race of the BH (solid red). Comparing 
Mc to the BH accretion rate M^isc demonstrates that during 
the peak episodes of accretion Mviac ~ Mc, so that radia- 
tion pressure becomes dynamically important. Although it 
is certainly possible to have accretion rates smaller than Mc 
when there is insufficient gas to fuel the AGN, the accretion 
rate is limited to a maximum value of ~ Mc. 

Fig. [2] shows the surface density of gas within Race = 
4e = 0.19 kpc for the fiducial simulation and for a higher 
gas fraction simulation with fg — 0.3. As implied by Fig. [TJ 
there are two main epochs during which significant gas is 
driven into the nuclei of the galaxies: after first passage and 
at final coalescence. The physical origin of these high nuclear 
gas densities are, however, somewhat different. 

iMihos fc Hernguistl l| 19961 ) showed that the presence of 
a bulge like that in our simulation suppresses a nuclear star- 
burst after first passage during galaxy mergers, because the 
bulge inhibits the non-axisymmetric modes that drive in- 
fiow. In our fiducial simulation, the majority of the increase 
in star formation after first passage is due to gravitational in- 
stability and fragmentation of the gas, which produces dense 
regions of rapid star formation. Fig. [3] (left panel) shows 
the gas density in the vicinity of one of the incoming black 
holes at t = 0.74 Gyr, midway through the first peak in star 
formation; the companion galaxy is well outside of this im- 
age. Two knots of dense gas are clearly seen, both of which 
will soon enter Race, the BH accretion and feedback region. 
These two clumps are not the only ones that form after first 
passage, but they are the only clumps that survive to enter 
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Figure 3. Gas density in the vicinity of the BH for the fiducial simulation at t = 0.73 Gyr (left panel), just prior to the onset of 
significant BH accretion after the first close passage of the two galaxies, and t = f .71 Gyr {right panel), the peak of star formation and 
BH accretion after the galaxies and BHs have coalesced. The times of these images are labeled with blue circles in Figure [T] In the left 
panel, the image is for the less inclined galaxy and the companion galaxy is well outside the image. The images are 5.7 kpc on a side and 
brighter color indicates a higher density. The dark region in the center of each image is within Race of the BH and is evacuated by BH 
feedback. In the image just after first passage [lejt panel), the two bright white regions are gaseous/stellar clumps that fragmented by 
Toomre instability during first passage and then spiraled into the nucleus, fueling star formation and BH accretion. At final coalescence 
{right panel), the nuclear gas densities are significantly higher (see also Fig. [2]l and most of the gas resides in a ~ 1 kpc diameter d isk 
driven into the nucleus by non-axisymmetric stellar torques during the merger. These images were made using SPLASH l|Pricell2007l) . 
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Figure 2. The mean gas surface density S interior to the accre- 
tion radius Race = 4e = 0.19 kpc for the fiducial simulation with 
initial gas fraction /g = 0.1 (solid) and for the simulation with 
fg = 0.3 (dashed; run fg). 

the central region surrounding the BhQ Fig. [3] (right panel) 
also shows an image of the gas density in the nuclear re- 
gion &t t = 1.71 Gyr, near the peak of star formation and 
BH accretion and after the galaxies and BHs have coalesced. 
At this time, the gas density in the nuclear region is signifi- 

^ In the simulation with a higher initial gas density {fg = 0.3), so 
many fragments form at large radii and spiral into Race that the 
surface density in the central region remains elevated from first 
passage until the merger completes at t ~ 1.8 Gyr (see Fig. [2]|. 



cantly higher than at first passage (see also Fig. and most 
of the gas resides in a ~ 1 kpc diameter disk. This nuclear 
gas concentration is the diffuse ISM driven in from larger 
radii by non-axisymmetric stell ar torques during the merger 
(e.g.. lMihos fc Hernauis^ll996l ). 

The galaxies in our fiducial simulation are stable when 
evolved in isolation. The merger itself drives the gas to 
fragment by locally exceeding the Jeans/Toomre mass. In 
reality, the gas in such clumps might disperse after ~ a 
Myr be cause of stella r feedback not included in our calcu- 
lations (|Murrav et al.il2010i ). This would probably not sig- 
nificantly change our estimate of the star formation rate 
since we are already normalized to the observed Kennicutt 
relation; however, such dispersal would lead to little inflow 
of gas associated with the inspiral of stellar clusters and 
thus would suppress the flrst peak in BH accretion (see 
iHopk ms fc Quataertll2009l for a more detailed discussion). 
In ij3.3l we will return to these issues and show that the 
total stellar mass and BH mass formed during the merger 
are relatively insensitive to the details of our assumed ISM 
model. 

Fig. 2] shows the surface density of gas in the fiducial 
simulation (top panel) and for the run without feedback 
(bottom panel) as a function of distance from the BH at 
four times; the initial condition {t = 0), shortly after the 
first close passage of the two galaxies (t = 0.85 Gyr), near 
the peak of accretion {t = 1.71 Gyr) and at the end of the 
simulation {t — 2.85 Gyr). Once M ~ Mc at first passage 
~ 0.85 Gyr, gas is driven out of the nuclear region by the 
AGN's radiation pressure. Since at the same time gravita- 
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Figure 4. Comparison of gas surface density (= Mg[< r]/iTr^) 
versus distance from the BH in the fiducial simulation with feed- 
back (top) and without feedback (bottom). Four times are shown: 
t = 0, 0.85 Gyr (first passage), 1.71 Gyr (peak accretion), and 
2.85 Gyr (end of simulation). Note that the gas tends to pile up 
at Race = 0.190 kpc (shown by the vertical line) in the top panel. 



tional torques continue to drive gas inwards, the gas begins 
to pile up at ~ Race- The particular radius at which the 
pile up occurs of course depends on our choice of Race, and 
so the particular size of the evacuated region should not 
be taken too seriously. Qualitatively, however, the behavior 
in Fig. |4] is reasonable: the AGN pushes on the gas in its 
neighborhood until it deprives itself of fuel. 

Near the peak of activity at t = 1.71 Gyr, the gas sur- 
face density in the central Race — 0.19 kpc is a factor of 
~ 10 — 30 larger in the simulations without feedback (bot- 
tom panel of Fig.|3|. However, the gas density at large radii 
~ 0.5 kpc is not that different. The radiation pressure force 
from the BH thus largely affects gas in its immediate envi- 
ronment, rather than the entire gas reservoir of the galaxy. 
Another indication of this is that the star formation rate is 
very similar in the simulations with and without feedback 
(middle panel of Fig. [T]). 

3.2 Dependence on Parameters of the BH Model 

The models for BH accretion and feedback used here contain 
uncertain parameters. We have defined the three relevant 
parameters a, r, and Race in i]2.4l and motivated our fiducial 
values, but it is important to explore how our results change 
with variations about our fiducial parameters. 

The value of a parameterizes the efficiency with which 
gas accretes from ~ Race ~ 190 pc to smaller radii, en- 
capsulating both the efficiency of angular momentum trans- 
port and the effects of star formation on unresolved scales. 
Naively, a higher value of a would lead to a more mas- 



sive BH. This is, however, not the case, because during the 
epochs when the BH gains most of its mass, the accretion 
rate is set by the efficiency of feedback (eq. [H)) not by the 
available mass supply (see Figs[T]&|4|. To demonstrate this 
more explicitly, the top left panel of Fig. [5] compares the 
BH accretion rates for three simulations with feedback, but 
differing values of a (0.05, 0.15, and 0.3), to the simula- 
tion with no feedback, which has a = 0.15. The accretion 
histories for the three values of a are nearly identical. By 
contrast, the accretion rate is in general much larger in sim- 
ulations that neglect feedback (and is oc a). In addition to 
the constant a runs, we tested a model in which a was time 
variable, set by the local gas fraction near the BH (fidafg2 in 
Table [l}: a = 3/g , with fg determined within Race (in prac- 
tice a varied from ~ 2 x 10~* — 0.3). Although this precise 
functional form is somewhat arbitrary, such a variation is 
motivated by analytic arguments and numerical simulations 
which show that instabilities due to self-gravity do minate 
the t r ansport o f gas from ^ 100 pc inward ( Shlosma n et al.l 
Il990l : fifopkins fc Quataertll2009| ). For our a = 3fg simu- 
lation, we find that the peak accretion rates and final BH 
mass are very similar to the constant a simulations. This 
is consistent with our conclusion that in the limit of large 
fuel supply, feedback, rather than the efficiency of angular 
momentum transport, sets the rate at which the BH grows. 

The parameter r describes the efficacy of the feedback 
for a given AGN luminosity. The bottom left panel of Fig. 
[5] compares the BH accretion rate for the fiducial run with 
r = 10 (black) and a simulation with a smaller value of r = 3 
(orange). To the extent that the accretion rate is feedback 
limited and set by Mc in equation [8] M should decrease 
with increasing r. Physically, this is because larger r leads 
to a larger feedback force, which then requires a smaller 
accretion rate to provide the luminosity necessary to drive 
away the surrounding gas. This expectation is borne out by 
the simulations. To compare the numerical results with the 
scaling in equation [HI the bottom left panel of Fig. [5] also 
shows M for the fiducial simulation scaled by a factor of 
10/3 (dashed line). This scaled M of the fiducial simulation 
is in reasonably good agreement with the t — 3 simulation, 
particularly at the first and second peaks in M, when most 
of the BHs mass is accumulated. This demonstrates that the 
value of T does not significantly affect any of the qualitative 
behavior of how the BH grows, although it does determine 
the overall value of the BH mass. 

In the majority of the simulations presented here, the 
size of the region over which we apply the feedback and 
average the gas properties to calculate M, Race, is set to 
4e. The rationale for this choice was given in ^2Al but it is 
important to consider the effects of changing this value. The 
top right panel of Fig.[S]shows the mass accretion rate for the 
fiducial simulation and a simulation with Race = 8e — 380 
pc. The peak values of M and the time of the first and second 
peaks are reasonably similar in the two cases. The principle 
difference is that in the simulations with the larger value 
of Race, the feedback is less effective at clearing gas out of 
the nuclear region (because the force is distributed over a 
larger number of particles); this allows a higher level of A/ to 
be maintained after the first passage and final coalescence. 
We suspect that the fiducial simulation better approximates 
what a higher resolution calculation with radiative transfer 
would find, but this remains to be demonstrated. 
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Figure 5. Top Left: Comparison of the accretion rates for the run without feedback (red; fidNof, a = 0.15), and three runs with feedback: 
the fiducial simulation with a = 0.05 (black), the run with a = 0.15 (green; fidSa) and the run with a = 0.3 (blue; fid6a). Top Right: 
The accretion rate for the fiducial run (black) and the run with Race = 8e (magenta; fidSeps). Bottom Left: The accretion rate for the 
fiducial run (black) and the run with r = 3 (orange; fidTau). Also shown is M^isc for the fiducial run increased by a factor of 10/3 
(dashed line), as expected from eq. Bottom Right: The integrated black hole masses for all the runs in this Figure with t = 10. 



1000 



_ 100 r 

T - 

s 10 . 

1 r 




0.1 I • 1 • 1 • 1 • 1 • 1 • 1 

0.5 1 1.5 2 2.5 3 

t [Gyr] 

Figure 6. The star formation rate for the run with no feedback 
(red) and for runs with various values of the BH accretion and 
feedback parameters: a = 0.05,0.15,0.3 (black, green, blue), a = 
3/g (grey), r = 3 (orange), and Race = 8e (magenta). All of these 
models have very similar star formation histories. 



The bottom right panel of Fig. [5] shows the integrated 
BH mass as a function of time for the fiducial simulation 
and for the variations in the feedback/accretion model con- 
sidered in this subsection that have the same value of r (but 
different values of a and/or Race)- The key result is that 
in the presence of feedback (all but the top curve), there is 
only a factor of ~ 3 change in the BH mass due to differ- 



ences in how we treat BH accretion and feedback. A factor 
of 6 change in a leads to only a 42% change in the final BH 
mass. This is because most of the BH mass is gained during 
the final coalescence of the two galaxies, at which point the 
BH accretion self-regulates and reaches the Eddington-like 
value in equation ([8]). The run without feedback (top curve), 
by contrast, has a factor of ~ 10 larger BH mass and the 
BH mass would scale linearly with the assumed value of a. 

The star formation rates for the simulations with differ- 
ent BH feedback parameters are all shown in Fig. [6] (this in- 
cludes the fiducial simulation with and without feedback and 
the runs with a = 0.15, 0.3, 3/g , t = 3, and Race ~ 8e). This 
figure demonstrates that the precise parameters of the BH 
feedback model have little effect on the galaxy-wide proper- 
ties such as the star formation rate: the total mass of stars 
formed in simulations with different BH feedback parame- 
ters differ by less than 5%. 

In previous simulations of BH growth and feedback, 
AGN feedback acting on dense gas in galaxies has be en in- 
voked to quench star formation (jSpringel et al.l l2005aV Our 
results demonstrate, however, that this is by no means guar- 
anteed (we refer here to 'quasar' feedback on cold dense gas, 
not the effect of AGN on hot dilute gas in galaxy groups and 
clusters). In our calculations BH growth is self-regulated and 
closely connected to the properties of the surrounding galaxy 
(e.g., eq. [8|. However, the BHs dynamical influence is cen- 
tered in the galactic nucleus (< 300 pc); as a result, the BH 
does not significantly alter the star formation history during 
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a merger. In this scenario, the merger remnant can nonethe- 
less be relatively quiescent ("red and dead") because the 
burst of star formation uses up much of the available gas. 

3.3 Effects of the ISM Model 
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Motivated by observations (e.g.. [Pownes fc Solomoiilll998l l, 
we have reduced the effective sound speed in GADGET'S 
subgrid ISM model (see i ]2.2p . There is nonetheless consid- 
erable uncertainty in the accuracy of this (or any other) sub- 
grid model. To study in more detail the effects of the ISM 
model on our results, we performed two additional simula- 
tions at our fiducial galaxy mass with the subgrid interpola- 
tion parameter Qeos = 0.2 and 0.07 (see eq. 3)), and without 
the factor of 10 reduction in pressure used in our fiducial 
simulation (an additional simulation with Qcos ~ 0.07 at a 
lower galaxy mass will be discussed in ijjllj The three dif- 
ferent ISM models have Cs and Q within a factor of ~ 2 of 
one another at all radii, with the Qbos ~ 0.2 model having 
the largest values of Cs and Q, and our fiducial model having 
the smallest values. The parameter Q is initially ~ 3 for our 
fiducial simulation at the disk scale length Ra, which is why 
the merger can induce significant fragmentation of the gas 
(Fig. [3|. Given the limited physics included in the subgrid 
model, we do not believe that it is feasible to unambigu- 
ously conclude which of these ISM models is more realistic. 
These models thus provide an indication of the systematic 
uncertainty introduced by our treatment of the ISM. 

Fig. [7] compares the BH accretion history (top panel), 
the star formation rate (middle), and the integrated BH 
mass and mass of new stars formed during the merger (bot- 
tom) for the three runs with differing ISM models. For both 
the fiducial run and the run with qeos = 0.07 there is sig- 
nificant fragmentation after first passage, which generates 
the first peak in star formation and BH accretion. By con- 
trast, the run with qsos = 0.2 shows no evidence for gas 
fragmentation or a pronounced peak in activity at first pas- 
sage. Despite these differing initial histories, the final result 
of the merger is very similar in all three cases: the large star 
formation rates and BH accretion rates coincident with the 
final coalescence of the two galaxies are not due to fragmen- 
tation, but are instead largely due to the inflow of diffuse gas 
to smaller radii. Moreover, the final BH mass and the total 
amount of new stars formed during the merger are similar in 
all three cases. Thus, despite uncertainties in the model of 
the ISM, we find relatively r obust integrated quantiti es (as 
did the earlier calculations o f lHernauist fc Mihoslll995l ). The 
precise time dependence of the star formation and BH ac- 
cretion (i.e., the lightcurves) are, however, significantly more 
uncertain and sensitive to the details of the model. 

3.4 Galaxy Parameters 

Having shown that the final BH mass and new stellar mass 
do not depend strongly on the uncertain parameters in our 
accretion, feedback and ISM models, we now examine how 
our results vary with galaxy properties such as the total 
mass, gas fraction, and bulge-to-disk ratio. 



3 We used Tsn = 4 x 10* K, Aq = 4000 and t° = 8.4 Gyr for 
these calculations; these values are different from those in our 
fiducial simulation, and are chosen to fix the total star formation 
rate for our isolated fiducial galaxy at IA/q yr"-'^. 
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Figure 7. Comparison of three simulations that differ only in 
the ISM models: fiducial (black), Qeos = 0-2 (red), and Qeos = 
0.07 (blue). The panels show the viscous accretion rate (top), 
star formation rate (middle), and the integrated black hole mass 
and mass of new stars formed (bottom). The three different ISM 
models have Cs and Toomre Q within a factor of ~ 2 of one 
another at all radii; the Qcob = 0.2 model has the largest values 
of Cs and Q and our fiducial model has the smallest values. 



Fig.[8]shows the BH accretion histories (top panel), star 
formation rate (middle), and integrated BH mass (bottom) 
for four runs with different total galaxy mass. The mod- 
els cover a factor of 30 in galaxy mass, from 0.1-3 times 
our fiducial mass. The BH and star formation parameters 
are identical in the four simulations, while the gravitational 
force softening and structural parameters (e.g., disk scale 
length, bulge radius) change with the total mass (see i]2.ip . 

Fig.[8]shows that the BH accretion rates and integrated 
BH masses increase with galaxy mass as expected from equa- 
tion[Sl However, there is a clear difference between the lower 
and higher mass simulations: the two higher mass simu- 
lations show evidence for the first peak in star formation 
and BH growth that we have shown is due to fragmentation 
near first passage, while the lower mass runs do not. This is 
largely a con sequence of the fact that observed disks have 
Rd oc M^/^ ()Shen et al.ll20oi l. so that more massive galax- 
ies have higher surface densities and are thus more suscepti- 
ble to gravitational instability (our ISM model counteracts 
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Figure 8. Comparison of four simulations that differ only in 
the total galaxy mass: Mfi^ = 3.88 X lO^^^Mo (fiducial; black), 
^Mfia (red), O.^Mfia (blue), and O.lMjid (magenta). The three 
panels show the viscous accretion rate (top), star formation rate 
(middle), and the integrated black hole mass (bottom). The same 
parameters are used in the BH accretion and feedback models. 



this slightly, but not enough to stabilize the higher mass 
disks). It is important to reiterate, however, that modest 
changes to the subgrid sound speed can change whether or 
not the gas fragments near first passage ( H3.3P so it is not 
clear if the difference as a function of mass in Fig.|8]is robust. 

In addition to the systematic change in the importance 
of fragmentation near first passage. Fig. [S] also shows dif- 
ferences in the late-time BH accretion between the low and 
high mass simulations. In particular the two smaller mass 
runs each show a period of increased accretion after the 
main peak during the merger. In these cases the new stars 
formed around final coalescence develop a bar in the inner 
~ Race of the galaxy. This helps drive some of the remain- 
ing gas into the accretion region leading to the increased 
accretion at late times. There is a milder version of this late- 
time accretion in the fiducial mass geos = 0.2 model without 
fragmentation in Fig. [T] Interestingly, there is no analogous 
late-time inflow of gas to within Race in our low mass galaxy 
simulations without BH feedback. The late-time activity is 
also particularly sensitive to the accretion model at a time 
when the non-axisymmetry produced by the merger has died 



away (so that a may in reality decrease significantly). For 
these reasons, we regard the late time growth in Fig. [8] as 
an interesting deviation from self-similarity in the dynamics, 
but not necessarily a particularly robust one. One important 
point that this highlights, however, is that because our im- 
plementation of BH growth and feedback does not unbind a 
significa nt amount of cold gas at late times (unlike calcula- 
tions bv lSpringel et al.ll2005al ). the predictions of our model 
are more sensitive to the post galaxy coalescence physics. 

In addition to the fiducial gas fraction [fg — 0.1) simu- 
lations that we have largely focused on, we performed simu- 
lations with an initial gas fraction of /g = 0.3 for our fiducial 
galaxy mass and at one tenth this mass. The qualitative dif- 
ference in behavior with galaxy mass in Fig. |8] persists in 
the higher gas fraction runs. In particular, in the low mass 
/g — 0.3 simulation, the gas does not fragment, while it does 
in the higher mass fg = 0.3 simulation. Fig. [2]- discussed in 
>i3.1l - explicitly shows the increase in the gas surface density 
within Race produced by this at early times. 

A final property of the galaxy model that we varied was 
the bulge to disk mass ratio. The majority of our runs in- 
clude a bulge with one third the mass of the disk; we also 
ran one simulation with an initial bulge of one fifth the disk 
mass, at the fiducial galaxy mass. The final BH mass and 
total mass of stars formed differ by less than 3% each com- 
pared to the fiducial simulation. 

4 THE Mbh - o CORRELATION 

Previous numerical studies using models of BH growth and 
feedback different from those considered here have repro- 
duced a number of the observed cor relations between mas- 
sive BHs and their host galaxies (e.g. , iDi Matteo et al.]|2005l; 
ISazonov et"aLl l2005l : I Younger etal] |2008| ). I Younger et"al] 
( 20081 ) argue that the galaxy-BH correlations in simulations 
(in particular, the BH fundamental plane) are relatively in- 
dependent of the trigger of BH growth, be it minor mergers, 
major mergers, or global instabilities of galactic disks. Based 
on the calculations to date, however, it is unclear to what 
extent the simulated BH-galaxy correlations depend on the 
details of the BH feedback or accretion models. In this sec- 
tion we assess this question by quantifying the Mbh - cr 
relation produced in our models. 

We define o of our model galaxies using a method analo- 
gous to that of observers: we first project the mass density of 
the stellar particles into cylindrical bins, and compute the 
half-mass(light) radius Re- We then compute the velocity 
dispersion weighted by the surface brightness via 



JR„, 



„2 



{R)I{R)RdR 



I{R)RdR 



(9) 



where I{R) is the projected 2-d stellar mass profile, aios is 
the line of sight velocity dispersion, and Rmin ~ 2e to ensure 
that there are that no artificial effects introduced by the 
force softening. We repeat this calculation along 1000 lines 
of sight with random viewing angles through the center of 
mass of the merger remnant. The a quoted in this paper and 
listed in Table 1 is the median value over the 1000 lines of 
sight. 

Fig- m shows the correlation between the final BH mass 
Mbhj and the a of the merged galaxy for most of the simu- 
lations in Table [1] different total galaxy masses (black), dif- 
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Figure 9. The Mgu j-a relation for the simulations in this pa- 
per, along wit h the observed r e lation (solid) and one sigma scatter 
(dotted) from lGiiltekin et alj |20oi). The final BH mass Mbhj 
in all of the simulations has been linearly scaled to r = 25 from 
the value used in the simulation (typically r = 10). The simula- 
tions are generally quite consistent with observations; we do find 
indications of a slight flattening in hlBH.f — cr at low BH masses. 



ferent values of the accretion parameter a (red circle), alter- 
nate ISM models (red x), higher gas fraction (blue square), 
alternate bulge mass (red square), difTerent values of r (blue 
circle), and the resolution studies in Appendix A (grey). The 
solid line indicates the mean relation from the compilation of 
observational results in iGiiltekin et al] (| 2009') while the dot- 
ted lines are the 1 — a error bars. We have linearly rescaled 
all of our final BH masses to a value of r = 25, using the 
fact that both the analytic and numerical results are con- 
sistent with Mvisc and Mbhj being oc r~^. The value of 
r = 25 is chosen so that the rescaled fiducial simulation lies 
appro ximately on the Mbh — cr relation of Giiltekin et ahl 
(|2009l ). For our fiducial simulation carried out with r = 3 
and T = 10, a linear scaling of Mbhj with is accurate 
to about 2% (e.g., Table [T] and Fig. [S]). We also carried out 
our fiducial simulation with r — 25; this is consistent with a 
linear scaling of Mbhj from r = 3 to ~ 50% (Table [T]). For 
nearly all of our simulations, rescaling to r = 25 amounts to 
dividing the final BH mass by a factor of 2.5. 

Previous analytic arguments were able to reproduce the 
Mbh — (J relatio n with t ~ 1, rather than requiri ng t ~ 25 as 
we do here (e.g.. lKinBir2003l ; iMurrav et al.ll2005l '). These cal- 
culations, however, assumed /g = 0.1. While perhaps appro- 
priate on average, this is not appropriate in galactic nuclei 
where the gas densities are higher. The analytic derivations 
also assumed that BH growth terminated when the system 
reached the luminosity (accretion rate) at which radiation 
pressure balances gravity (eq. [8|. In reality, however, the 
luminosity must exceed this critical value by a factor of sev- 
eral in order for gas to be efficiently pushed around in the 
galactic nucleus (as shown explicitly in the test problems in 
the Appendix). Moreover, the BH continues to accrete some 
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Figure 10. The ratio of the final BH mass to the BH mass at 
the peak of accretion for the simulations in Fig.[9l using the same 
symbol types. This quantifies the extent to which late-time accre- 
tion increases the BH mass. The late-time increase in BH mass 
for many of the lower mass systems produces the slight flattening 
in Mbh — c in Fig. [9]at low masses; see the text for a discussion 
of the robustness of this result. 



mass even after reaching M^- Fig. [TO] shows this explicitly 
via the ratio of the final BH mass to the BH mass at the 
peak of activity for all of the simulations in Fig. l9p| The 
net efi'ect of the differences between our simulations and the 
simple analytic calculations is that a much larger feedback 
force per unit BH mass (r ~ 25, not ~ 1) is required for 
consistency with the observed Mbh — o" relation. The phys- 
ical implications of this larger value of t for models of AGN 
feedback will be discussed in § [5] 

The scatter in BH mass in Fig. [5] at our fiducial mass 
scale of (T ~ 175 km s~^ is reasonably consistent with the ob- 
served scatter. In the simulations we have varied the BH ac- 
cretion model (a), the ISM model, numerical resolution, size 
of the feedback/accretion region Race, and galaxy properties 
such as the total mass, gas fraction, and bulge to disk ratio. 
It is encouraging that all of these simulations produce BH 
masses within a factor of few of each other. The largest BH 
mass at CT ~ 175 km s"'^ is the simulation with an initial gas 
fraction of /g = 0.3; since this run has a larger gas density at 
small radii close to the BH (Fig. [2]), it should probably also 
have a larger r, which would reduce the BH mass further, 
in better agreement with the data. It is difficult to make 
this comparison to the observed scatter more quantitative 
given the limitation that our simulations are all equal-mass 
non-cosmological binary mergers on the same orbit. 

The numerical results in Fig. |9] suggest a slight flatten- 
ing of the Mbh — o" relation at cr < 100kms~^. This is in 
large part a consequence of the additional mass gained by 
the lower mass BHs after their peak of activity (see Figs. [S] 
& 1101 in particular the fiducial simulations labeled by black 
squares in Fig. llOp . This change in behavior at lower masses 
is primarily due to the fact that the lower mass galaxies are 
less prone to fragmentation than the more massive galaxies 
( >i3.4|l . Without the fragmentation after first passage, more 



To account for the fluctuating nature of the BH accretion rate 
in some of the simulations, we deflne the BH mass at "peak" to 
be the mass when M drops by a factor of 10 from its peak value. 
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gas is available to feed the BH at late times leading to the 
slightly higher BH mass. As discussed in i]3.4l it is unclear 
how robust this late time accretion is. In fact, a low mass 
galaxy simulation with an alternate ISM model (gcos = 0.07) 
does not show significant late-time accretion, leading to a 
BH mass in good agreement with the extrapolation from 
higher cr (red x at low mass in Figs. [9l fc I10|) . We thus re- 
gard the case for flattening of Mbh — o" at low masses in our 
models as somewhat tentative; our results may instead indi- 
cate enhanced scatter at low masses rather than a change in 
the mean relation. More comprehensive numerical studies of 
these lower mass systems will be needed to distinguish these 
two possibilities. 

5 DISCUSSION AND CONCLUSIONS 

We have presented a new method for simulating the growth 
of massive BHs in galaxies and the impact of AGN activ- 
ity on gas in its h ost galaxy (see also our related Letter; 
iDeBuhr et al]|2009D . In this method, we use a local viscous 
estimate to determine the accretion rate onto a BH given 
conditions in the surrounding galaxy (eq. [S]) , and we model 
the effect of BH feedback on ambient gas by depositing mo- 
mentum radially away from the BH into the surrounding gas 
(eq.EJ. 

Our accretion model qualitatively takes into account 
the angular momentum redistribution required for accre- 
tion of cold gas in galaxies and is thus more appropriate 
than the spherical accretion estimate that has been used 
extensively in the literature. In our feedback model, the ap- 
plied force is given by tL/c, where the AGN's luminosity 
L is determined by our BH accretion model, and the net 
efficiency of the feedback is determined by the total op- 
tical depth r of the galactic gas to the AGN's radiation, 
which is a free parameter of our model. Previous calcula- 
tions have demonstrated that only when the gas fraction 
in a galaxy decreases to < 0.01 can the A GN's radiation 
Com pton heat matter to high temperatures (|Sazonov et al.l 
l2005h . More generally, the cooling times in gas-rich galaxies 
are so short that the primary dynamical impact of the AGN 
on surrounding gas is via the momentum imparted by the 
AGN's outflows or radiation. It is thus not physically well- 
motivated to model AGN feedback by depositing energy, 
but not momentum, i nto surrounding gas, as many calcula- 



tions have done (e.g., Di Matteo et al. 20051: Springel et al.l 
' lOstriker et al.i (|2010l ) 



l2005al : iKawata fc GibsonI l2005h 



for related points. 

Throughout this paper, we have focused on BH growth 
d uring major m ergers of spiral galaxies. As demonstrated 
m iDeBuhr et al] (j2Q09), our model leads to a self-regulated 
mode of BH accretion in which the BH accretion rate is 
relatively independent of the details of the BH accretion 
model (see Fig. [SJ . This is because the accretion rate self- 
adjusts so that the radiation pressure force is comparable to 
the inward gravitational force produced by the host galaxy 
(see eq. |5]|. This self-regulated mode of BH accretion is a 
robust feature of all of our simulations during periods of 
time when there is a significant nuclear gas reservoir - it 
thus applies precisely when the BH gains most of its mass. 

One important consequence of this self-regulated ac- 
cretion is that AGN feedback does not drive significant 
large-scale outfiows of gas (in contrast to the models of 



ISpringel et"alll2005al ). For example, the surface density pro- 
files in Fig. |4] show that AGN feedback causes gas to pile up 
at a few hundred pc rather than being completely unbound 
from the galaxy - this precise radius should not be taken 
too literally since it is a direct consequence of the fact that 
we implement feedback and determine the BH accretion rate 
only within a radius Race ^ few hundred pc. Nonetheless, 
we believe that this general result may well be generic: be- 
cause the BH accretion rate is determined by the gas content 
close to the BH, the AGN can shut off its own accretion be- 
fore depositing sufficient energy to unbind all of the gas in 
the galaxy. If we artificially hold the luminosity of the AGN 
constant in time at a value exceeding the critical value in 
equation Q, then the AGN does eventually unbind all of 
the surrounding gas (see, e.g.. Figs. IBll fc lB3l in Appendix 
B). However, both our isothermal sphere test problem (Fig. 
IB6|) and our full merger calculations show that when the 
BH accretion rate is self-consistently determined by the gas 
properties in the central 100 pc of the galaxy, the AGN 
simply never stays 'on' long enough to unbind all the gas. 

Our results do not, of course, preclude that AGN drive 
galactic winds. For example, some gas may be unbound by a 
high speed wind/jet produced by the central accretion disk 
(which is not in our simulations). In addition, at later stages 
of a merger or at large radii the gas fraction can be suffi- 
ciently low (< 0.01) that gas can be Gompton h eated to high 
temperatures and potentially unbound (e.g., ICiotti et al-l 
I2OIOI ). This may in fact be sufficient to quench star forma- 
tion at late times, but only once most of the gas has already 
been consumed into stars (so that fg < 0.01). Our results 
do suggest that AGN feedback does not quench star forma- 
tion by unbinding a significant fraction of the cold dense 
gas in a galaxies int erstellar medium (in contrast to, e.g., 
ISpringel et al.|[2005al ). In future work it will be important 
to assess whether variability in t he accretion rate on smaller 
scales than we can resolve (e.g., iHopkins fc Ouataertll2009l : 
iLevine et ai1l2010^ modifies this conclusion; such variability 
could produce some epochs during which AGN feedback has 
a significantly larger effect on the surrounding gas. Another 
improvement would be to carry out radiative transfer calcu- 
lations and assess what fraction of the AGN's radiation is 
absorbed at large radii in a galaxy (~ kpc) where the gas 
has a lower surface density and is thus easier to unbind. 

Our simulations cover a factor of ~ 30 in galaxy mass. 
The final BH mass in our calculations is oc r^^ since a larger 
value of r corresponds to a larger momentum deposition per 
unit BH mass. We find reasonable consistency with the nor- 
malization of the observed Mbh — rela tion for t ^ 25. 
To compare this result to previous work bv lDi Matteo et al.l 
(2005), we note that a momentum deposition of P corre- 
sponds to an energy deposition rate of i5 ~ Pa when the 
feedback is able to move the gas at a speed comparable to 
the velocity dispersion a (which is required for efficient self- 
regulation of the BH growth). For r ~ 25, our results thus 
correspond to £ ~ 25 L (j/c:^ 0.02 L (g/200 km s~^). This is 
similar to the results of iDi Matteo et al] l|2005h , who found 
that depositing ~ 5% of the BH accretion energy in the sur- 
rounding gas was required to explain the Mbh — o" relation. 
It is encouraging that these two different sets of simulations, 
with different BH accretion and feedback models, agree at 
the factor of ~ 2 — 3 level on the energetics required to 
reproduce the Mbh — cr relation. 
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The value of r ~ 25 required to explain the normal- 
ization of the Mbh — o" relation has strong implications 
for the dominant physics regulating BH growth. The sim- 
plest models of super-Eddington winds from an accretion 
disk close to the BH are ruled out because they typically 
have r ~ 1, i.e., a mo mentum flux comparable to that of the 
initial radiation field (|King| [ioOSl l . Similarly, the radiation 
pressure force produced solely by the scattering and absorp- 
tion of the AGN's U V radiation by dust corresponds to r ~ 1 
l|Murrav e t al. 2005) and is thus not sufficient to account for 
the level of feedback required here. Rather, our results sug- 
gest that most BH growth happens when the nuclear regions 
are optically thick to the re-radiated dust emission in the 
near and far-infrared, so that t ^ 1. This is consistent with 
observational evidence in favor of a connection between BH 
growth, quasars and luminous dust-enshrou ded starbursts 
such as ULIRGs and sub-mm galajcies (e . g., Sanders et al.l 
1 19881 : iDasvra et al]|2006l : [Alexander et al.l 120081 '). Quantita- 
tively, the observed stellar densities at radii 1 — 10 pc in 
elliptical galaxies reach ~ 20 g cm' -2 l|Hopkins et al.ll201(t ). 
implying r ~ 100 if a significant fraction of the stars were 
formed in a single gas-rich epoch. It is encouraging that 
this is within an order of magnitude of (and larger than!) 
the value of r we find is required to explain the observed 
Mbh — o" relation. 

A fixed value of r ~ 25 independent of galaxy mass 
produces an Mbh ~ o" relation with a slope and scatter in 
reasonable agreement with observations (see Fig. [9]). Assess- 
ing the scatter more quantitatively will require a wider sur- 
vey of merger orbits. We do find some tentative evidence 
for a shallower slope in the Mbh — o" relation at the lowest 
galaxy masses, corresponding to ct < lOOkms"^. This range 
of masses is precisely where the observational situation is 
particularly unclear, with, e.g., possible differences between 
the BH -galaxy correlation s in classical bulges and pseudo- 
bulges iGreene et aLllioO^ V It is also unclear whether major 
mergers are the dominant mechanism for BH g rowth in these 
lower mass galaxies (e.g.. lYounger et al.lb'oog '). 

Our simulations show that fragmentation of a galactic 
disk into clumps can be efficiently induced by a merger (e.g.. 
Fig- El, even when an isolated galaxy with same properties 
is Toomre stable (see, e.g., Wetzstein et al, 2007 ; for related 
ideas in the context of dwarf galaxy formation in tidal tails). 
As Figure [7] demonstrates, this fragmentation can produce 
a significant increase in star formation during the first close 
passage of galaxies even when there is little inflow of the 
diffuse ISM (because suc h inflow is suppressed by a bulge 
until later in the merger; iMihos fc HernauistllT99g ). In our 
simulations we often see a corresponding increase in the BH 
accretion rate due to the inspiral of dense gas-rich clumps 
(Fig- El . The inflow of gas by this process may, however, be 
overestimated: stellar feedback not included in our simula- 
tions can unbind the gas in star clusters on a timescale of 
^ a Myr, retu rning most of the gas to the diffuse ISM (e.g., 
iMurrav et al] 2010; Hopkins & Quatacrt 2009). 

Our calculations use subgrid sound speeds motivated by 
the observed turbulent velocities in galaxies ( H2.2p . We thus 
believe that our ISM model is physically well-motivated, 
even though the use of a subgrid sound speed necessarily in- 
troduces some uncertainty. Overall, the presence/absence of 
large-scale clumping of the ISM does not signiflcantly change 
the final BH mass or the mass of new stars formed in our 



simulations. It can, however, change the star formation rate 
and BH accretion rate as a function of time, particularly 
near the first close passage during a merger. 

The tentative change in the Mbh — o" relation we find 
for lower mass galaxies is largely due to our treatment of the 
ISM, rather than our BH feedback or accretion model. For 
a given gas fraction, lower mass galaxies have a lower gas 
surface density and thus the ISM is less prone to fragmenta- 
tion ( H3.4l and Fig. [5|. Without the fragmentation after first 
passage, more gas is available to feed the BH at late times 
leading to somewhat higher BH mass (Fig. llOp . 

The BH accretion and feedback models used in this pa- 
per can be signiflcantly improved in future work, allowing 
a more detailed compariso n to observations. For example, 
iHopk ms fc QuataertI (|2009l ) carried out a large number of 
simulations of gas infl ow in galactic nuc lei from ~ 100 pc 
to < 0.1 pc (see, e.g., iLevine et al.ll2010l for related work). 
These can be used to provide a more accurate estimate of the 
BH accretion rate given conditions at larger radii in a galaxy 
(Hopkins & Quataert, in prep). Another important improve- 
ment would be to use a radiative transfer calculation to self- 
consistently determine the infrared radiation fleld produced 
by a central AGN (and distributed star formation). This 
could then be used to calculate the radiation pressure force 
on surrounding gas, eliminating the need for our parameter- 
ization of the force in terms of the optical depth r. 
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APPENDIX A: RESOLUTION STUDIES 

In this section, we describe some of our resolution tests both 
with and without BH feedback. In the absence of feedback, 
the well-posed questions for resolution studies include both 
how the gas properties as a function of radius and time de- 
pend on the resolution and how integrated properties of the 
galaxy (e.g., the star formation rate) depend on resolution. 
However, the feedback, when present, has a nontrivial de- 
pendence on the resolution and it is by no means clear that 
the nonlinear system will in fact converge in a simple way 
with increasing resolution. Physically, e.g., the AGN's radi- 
ation pressure has the strongest effect on the gas that con- 
tributes the most to the optical depth, which is largely de- 
termined by the column density (the dust opacity being only 
a relatively weak function of temperature for the conditions 
of interest). Higher resolution simulations can resolve higher 
volume and column densities, largely at smaller radii close 
to the BH, and thus may change some of the details of the 
BH feedback. Indeed, Fig. [4] shows that the column density 
increases towards smaller radii in our simulations. 

We first consider the question of how the nuclear gas 
properties depend on numerical resolution in the absence of 
feedback. To this end, the top left panel of Fig. lAll shows 
the BH accretion rate Mvisc calculated for three different 
particle numbers Np = 1.6 x 10^4.8 x 10^ and 1.6 x 10*^, 
with the gravitational force softening e oc A'jJ^^Q To make 
a fair comparison, the accretion rate is evaluated within a 
fixed volume {R = 406 pc) and for a = 0.05 for all of the 
simulations. This choice corresponds to _R = 4e for the low- 
est resolution run, but is i? ~ 8.6e for our fiducial resolution 
simulation. Fig. lAll shows that the lowest resolution simula- 
tion (red) does not adequately resolve the fragmentation of 
the gas, and the resulting peak in the accretion rate, near 
first passage. The medium and higher (= our fiducial) res- 
olution simulations, however, agree reasonably well, except 

^ In Fig. lAll Mvisc for the simulations without feedback (upper 
left) is calculated from the simulation snapshots and the accretion 
rate is not Eddington limited. The data outputs were relatively 
infrequent and attempting to integrate the BH mass over such 
large timesteps was inaccurate. 



© 0000 RAS, MNRAS 000, 000-000 



BH Growth with Feedback by Radiation Pressure 17 




t [Gyr] t [Gyr] 



Figure Al. Top left: BH accretion rate for simulations without feedback at three resolutions (LRfidNof, MRfidNof, fidNof in red, blue 
and black respectively); M^isc was computed for a = 0.05 and using the same value of Race = 406 pc at all three resolutions (this 
corresponds to 4e for the lowest resolution). Top right: BH accretion rate with feedback at the same three resolutions using Race = 406 
pc. Bottom left: BH accretion rate with feedback at the same three resolutions using Race = 4e; here the accretion rate and feedback are 
calculated in different physical volumes at different resolutions. Bottom right: BH masses for all of the runs in this Figure with feedback. 
Solid lines are for Race = 4e while dashed lines are for Race = 406 pc. 



for a slight difference in the slope of Mviseit) at late times. 
Computed over a larger volume (~ kpc), the agreement be- 
tween these runs improves. 

To assess the convergence in the presence of feedback, 
the top right panel of Fig. lAll shows the BH accretion rate 
Mvisc evaluated just as in the top left panel, i.e., using a fixed 
Race = 406 pc, in simulations with the same three particle 
numbers and force softening. Again the lowest resolution 
(red) simulation is clearly not adequate, but the medium 
(blue) and high (black) resolution simulations agree well; 
the integrated BH mass differs only by 2% in the latter two 
simulations. 

As a final resolution test, the bottom left panel of 
Fig. lAll shows the BH accretion rate as a function of time 
in simulations with the same three resolutions and force 
softening, but in which Race = 4e. Thus in this case the 
accretion rate is determined, and the feedback applied, on 
increasingly small spatial scales in the higher resolution sim- 
ulations. This is probably the most physically realistic (see 
i|2.4|) . This panel shows that the large peak of accretion at 
final coalescence (t ~ 1.8 Gyr) is quite similar in all three 
cases. This is set by the physics of feedback by momentum 
deposition and is a robust property of all of our simulations. 
A corollary of this is that the final BH mass, as shown in 
the bottom right panel of Fig. lAll is the same to within a 
factor of ~ 2 for the three different resolutions. However, the 
results in the lower left panel of Fig. lAll also clearly demon- 



strate that the detailed evolution of the accretion rate is sen- 
sitive to the resolution. This is not particularly surprising: 
at fixed resolution. Fig. [5] has already demonstrated that the 
details of AIvise{t) depend on the value of Race - although, 
again, neither the integrated BH mass or star formation rate 
do. One implication of these results is that it is difficult for 
current simulations of BH growth to make quantitative pre- 
dictions about the Ught curves of AGN activity triggered by 
mergers. 



APPENDIX B: CODE VERIFICATION 

We have tested our modifications to GADGET on a num- 
ber of simplified problems that have answers that can be 
easily obtained through other methods. §B1 describes tests 
of the additional momentum feedback force applied to a thin 
spherical shell of gas. §B2 describes tests in which the force 
is applied to the gas particles in the central regions of an 
isothermal sphere. Two ways of implementing the force are 
tested: to a fixed number of particles around the BH, and 
to all particles within a fixed region Race around the BH. 

As we are concerned with the performance of our BH 
accretion and feedback model, in all of the tests presented in 
this appendix, t he multiphase equation of stat e and star for- 
mation model of ISpringel fc Hernguistl l|2003l ) are not used; 
instead we use an adiabatic equation of state with 7 — 5/3. 
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Figure Bl. Time evolution of the radius of the test shell for three 
values of radiation force: A = 0.5,1.0,2.0 (dashed curves). The 
results match closely with the solutions from integrating eq. IIB2II 
(superposed grey curves). Here the force is applied to the 25000 
innermost gas particles of the 5 X 10* that make up the shell. 

Time is in units of to = ^ Tq/GMbh and the radius is in units 
of ro, where ro is the initial radius of the gas shell. 
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Figure B2. Time evolution of the radius of the test shell for 
three values of N/N^hell - 0.5 (solid), 0.25 (dashed), and 0.1 (dot- 
dashed). The numerical solutions are normalized by the exact 
solution from eq. |)B2|| . The radiation force is fixed to be A = 
2.0. The radius r{t) changes by only about 1% as N is changed, 
indicating that our results are insensitive to the exact number of 
particles to which the radiation force is applied. 



Bl Gas shells 

To test that the code is applying the radiation pressure force 
in equation (O correctly, we have run the code for a simple 
system containing a black hole particle with a large mass 
and a thin spherical shell of gas with negligible temperature, 
pressure and mass. As this gas resides in a thin shell, this 
problem is more well-posed if we apply the radiation force 
to a fixed number. A*', of gas particles. 

This system has a critical luminosity defined by the 
point at which the radiation force balances the inward pull 
of gravity. As the gas shell is of low temperature and pres- 
sure, we can ignore pressure forces. For a black hole of mass 
Mbh and a gas shell of mass m at a radius ro the critical 
luminosity Lc satisfies (we take t = 1 for simplicity) 



r ^MBHm 

Lc ~ (_r 5 C. 



(Bl) 



When the luminosity is set to this value, the gas shell should 
experience no net force. For other values of the luminosity, 
the expected behaviour can easily be calculated by noting 
that the gas shell, in the absence of any pressure forces, 
should have a radius, r(t), that satisfies 



cfrjt) 



GMsHm L 
r(t)2 



(B2) 



This is easily integrated to give the expected behavior of the 
gas shell. 

A number of tests of this system were performed with 
varying luminosities, parameterized by the ratio of the lu- 
minosity applied to the critical luminosity, 



(B3) 



Fig. IBll shows the exact result in grey, with the numerical 
solution from the modified version of GADGET in black, for 
runs with A — 0.5, 1.0 and 2.0. For these tests the number 
of particles in the shell is Nsheii = 50000, and the force was 
applied to A'^ = 25000 of them. In all cases, the numerical 



solution appears indistinguishable from the exact solution 
of eq. (|B2)) . 

We have also tested the dependence of the results on 
the value of N/Nsheii, the fraction of particles that receives 
the radiation force. Fig. IB2l shows the ratio of the numerical 
solution from our code to the exact solution for N/Ngheii = 
0.5,0.25, and 0.1 for the A = 2.0 model. This demonstrates 
that even though the magnitude of the force on an individual 
particle increases as N decreases, the overall dynamics of the 
shell is the same, with the radii differing by only ~ 1% in the 
three cases. This is primarily due to the fact that the SPH 
particles are coUisional and can thus transfer their motion 
to their neighbors via pressure forces. The extra momentum 
imparted to the subset of particles is transferred in part to 
the outer region of the shell, leading to the overall motion 
that agrees well with the exact solution. By extension, if N 
were to vary over the duration of the simulation, the results 
would also not depend strongly on the particular value. 



B2 Isothermal Sphere 

We have performed a second set of tests of the feedback 
model using an isothermal background given by a King 
model. The mass of the system is split into two parts. The 
bulk of the mass makes up the coUisionless background that 
is drawn from the full phase space distribution of the King 
model. A small fraction of the mass, fg — 0.05, is assigned 
as coUisional SPH particles. These gas particles follow the 
same spatial profile as the coUisionless background but are 
given zero initial velocities and a very low temperature. Both 
components are realized with 10^ particles. Finally, a black 
hole particle with a small mass is placed at rest at the center 
of the distribution. 

In the absence of feedback, the SPH particles are not 
in equilibrium by construction and should flow toward the 
center of the potential provided by the coUisionless back- 
ground. When the feedback is switched on in the isothermal 
King potential near the center, the feedback will again have 
a critical value set by force balance: 
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Figure B3. The density (top) and pressure profiles (bottom) in simulations of an isothermal sphere of gas embedded in an isothermal 
King potential. Three different simulations are shown: without feedback from a central black hole, A = 0.0 (left), A = 1.0 (middle) and 
A = 2.0 (right). Each simulation is shown at four times: the initial profile (black) and t = 0.16 (red), 0.32 (green), 0.48 (blue) Gyr. 



When the luminosity is below this value, we expect the ex- 
tra momentum to be insufficient to clear the gas out of the 
center. When the luminosity exceeds this value, the feed- 
back should be strong enough to clear the central regions 
of the distribution. To test this, we apply feedback with a 
constant luminosity. Again, we parameterize the strength of 
the feedback as X = L/Lc. 

We have tested two ways of assigning the radiation 
force. In the first case, the force is shared (equally) by a 
fixed number of gas particles nearest to the black hole. In 
the second case, the force is shared by all gas particles within 
a fixed radius of the black hole. We discuss the results sep- 
arately below. 

B2.1 Fixed N 

For the tests in this subsection, the radiation force is applied 
to a fixed number of gas particles: A'^ = 500. The King model 
has a = 100kms"\ */cr^ = 12 and a total mass of IO^^Mq. 

Fig. IB3I compare the density and pressure profiles of 
three runs with A = (i.e. no feedback; left panels), 1 (mid- 
dle), and 2 (right). Four timesteps are shown: t — (black), 
0.16 (red), 0.32 (green), and 0.48 Gyr (blue). As expected. 



the gas flows to the center in the absence of feedback, in- 
creasing the density and pressure as the gas begins to equi- 
librate in the background potential. The middle and right 
panels show that the feedback clearly has an effect on the 
gas at the center, providing some support for the incoming 
gas, allowing the gas to have a lower pressure. For the case 
with A = 2, the feedback is strong enough to effectively clear 
out the central region. 

The nature of the feedback allows a calculation of how 
the size of the evacuated region should grow with time. Ig- 
noring the thickness of the shell swept up as matter begins 
to be driven out by the feedback, momentum conservation 
gives 

-77 [MshM{r)dr dt] = ^ (B5) 

at c 

where Msheii{r) is the initial mass distribution of gas and 
Mbg{r) is the mass distribution of the background. Near the 
center of the initial distribution, both the gas and back- 
ground have an isothermal distribution, for which the mass 
increases linearly with the distance from the centre. This 
makes the right hand side of Eq. (|B5|l a constant. In this 
case, the size of the evacuated region, r{t), depends linearly 
on time: 

r{t) = ^2{\ ^ 1){1 - fg)at + C (B6) 
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Figure B4. Time evolution of the size of the evacuated region for 
the isothermal sphere test. The A = 4 simulation results shown 
(solid) match well with the analytic solution (Eg. IB6| dashed). 

where C is a constant of integration to account for the finite 
time required to form the sheU of swept up gas. 

Fig. lB4l shows the size of the evacuated region as a func- 
tion of simulation time for a run with A = 4, and the ex- 
act solution for a shell moving in an isothermal background 
(Eq. IB6|) with C chosen to match the position of the shell 
at t = 0.1 Gyr. The size of the evacuated region is defined 
by the position of the gas particle closest to the black hole. 
The agreement is very good with only slight deviation at 
the latest times. For the model employed, the potential is 
only isothermal near the origin, so when the shell expands 
sufficiently, the potential shallows and the shell should move 
faster than the prediction. This is indeed seen at late time 
in Fig. [Bl 

B2.2 Fixed Race 

For the galaxy merger simulations, we apply the force in- 
side a fixed Race throughout the simulation. In this section, 
we run a similar set of tests as in the previous subsection 
but we hold Race fixed. When the number of particles inside 
Race becomes small, however, the feedback force exerted on 
individual particles becomes spuriously large. We therefore 
impose an additional condition of minimum A'' on the feed- 
back. For the tests in this subsection, the feedback is applied 
to those particles inside Race, or to the innermost 100 gas 
particles if there are fewer than this inside Race- For the 
simulations in the main paper, however, there were always 
enough particles inside the accretion and feedback region to 
avoid the need for such a lower bound on A^. 

Our first test uses a constant L — 4Lc, and holds Race 
fixed. We use a King model as in the previous section, but 
with slightly different parameters to connect more closely 
to the our fiducial simulation: a = 160kms~^, "^/a^ = 12 
and a total mass of IO^^Mq. We tested this model for three 
different sizes of the accretion and feedback region: Race = 
0.7, 1.4 and 2.8 kpc. The smallest region has initially ~ 
500. Note that the values of Race used here are larger than 
those used in our galaxy merger simulations in the main 
text. These values of Race were necessary to ensure that 
Race contains a reasonable number of particles. In the galaxy 
merger calculations, the overall larger number of particles in 
the simulation and the high gas density in the central regions 
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Figure B5. Radius of the swept up shell for the isothermal sphere 
test with A = 4 and fixed Race- 0.7 (black), 1.4 (red), and 2.8 kpc 
(blue). To avoid numerical problems, the feedback was always 
applied to at least A'^ ~ 100 particles. The numerical results agree 
well with the dashed curve, which shows a numerical integration 
of the analytic equation for the shell radius (eg. IB5|I . 
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Figure B6. The black hole accretion rate for isothermal sphere 
simulations in which the full black hole accretion and feedback 
model are used (a = 0.1 and r = 1). Three different values of 
Race are shown: Race = 0.7 (black), 1.4 (red), and 2.8 kpc (blue). 
All three agree well with each other. 

imply that smaller values of Race can be reliably used. They 
are also more physical, as we argued in !j2]4l 

Fig. IB5I shows the position of the shell of swept up ma- 
terial for the three runs with Race = 0.7, 1.4 and 2.8 kpc 
in black, red and blue respectively. Initially, all the gas in- 
side Race experiences the extra force. As the region becomes 
more evacuated and the number of particles inside Race 
drops, we transition to applying the force to the A'^ = 100 
particles closest to the BH. The evolution of the shells in this 
case is quite similar to the evolution in the last section. The 
model used in this section is smaller in size and so the shell 
expands past the isothermal core of the King model earlier. 
As a result, it begins to accelerate outward sooner. How- 
ever, the tests with different Race have essentially identical 
evolution. 

Finally, we run a test in which we determine the lumi- 
nosity from the accretion rate as in Eq. (0, and increase 
the BH mass in time accordingly. This test thus employs 
the full feedback and accretion model of our galaxy merger 
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simulations. We use the same a = 160 km s^^ King model, 
and took a = 0.1 and t — 1 for the feedback parameters. 
The initial mass of the black hole was MBH,i = W^Mq. 

Fig. IB6I shows the accretion history of the BH for the 
runs with Race = 0.7, 1.4, and 2.8 kpc. In each test, the 
feedback is initially Eddington limited and it is not until 
about t = 0.3 Gyr that the luminosity approaches that re- 
quired to evacuate the gas out of Race- At this point, the 
gas begins to move out of Race and form a shell of material 
at 7? ~ Race- This shell then remains fairly steady as the 
accretion rate self-regulates around the critical luminosity. 
As the three values of Race are all inside the isothermal core 
of the King model, the critical luminosities (eq. IB4|) are the 
same, and we would thus expect the accretion rate to self- 
adjust to the same value at late times. This is indeed borne 
out in the simulations shown in Fig. IB6I Of these three runs, 
only the calculation with Race = 0.7 kpc spends a significant 
amount of time with fewer than 100 particles inside Race- 
Despite the large change in the size of the feedback region. 
Fig. IB6I shows that the evolution of the gas is quite similar. 
The black hole masses for these three runs differ by only a 
factor of ~ 2 at the end of the simulation. 
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